set.seed(15)
library(maftools)
library(survival)
library(survminer)

# 1. 数据加载 ----
laml.maf <- system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools') 
laml.clin <- system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools') 

# 2. 数据处理 ----
laml <- read.maf(maf = laml.maf, clinicalData = laml.clin)
samples <- laml@variants.per.sample$Tumor_Sample_Barcode
laml <- subsetMaf(laml, tsb = sample(samples, 150))

# 3. 绘制oncoplot ----
gene_summary <- getGeneSummary(laml)
top_genes <- c("TP53", gene_summary[1:14, Hugo_Symbol])

oncoplot(
  maf = laml, 
  top = 15, 
  genes = top_genes,
  fontSize = 0.8,
  titleFontSize = 1.2
)

# 4. 生存分析准备 ----
# 提取TP53突变样本
tp53_mutated <- laml@data[Hugo_Symbol == "TP53", unique(Tumor_Sample_Barcode)]
clin_data <- laml@clinical.data

# 处理生存时间（转换为年）
clin_data$time_years <- clin_data$days_to_last_followup / 365

# 修复生存时间中的无限值
valid_times <- clin_data$time_years[is.finite(clin_data$time_years)]
min_valid <- min(valid_times, na.rm = TRUE)

clin_data$time_clean <- ifelse(
  is.infinite(clin_data$time_years),
  min_valid,
  clin_data$time_years
)

# 5. 创建生存分析数据框 ----
surv_df <- data.frame(
  sample = clin_data$Tumor_Sample_Barcode,
  time = clin_data$time_clean,
  status = ifelse(clin_data$Overall_Survival_Status == "DECEASED", 1, 0),
  TP53 = ifelse(
    clin_data$Tumor_Sample_Barcode %in% tp53_mutated,
    "Mutated", 
    "Wildtype"  # 修复：添加野生型作为对照组
  )
)

# 6. 执行生存分析 ----
fit <- survfit(Surv(time, status) ~ TP53, data = surv_df)

# 7. 绘制生存曲线 ----
ggsurvplot(
  fit,
  data = surv_df,
  pval = TRUE,           # 添加p值显示
  conf.int = TRUE,
  risk.table = TRUE,
  title = "TP53 Mutation Survival Analysis",
  xlab = "Time (years)",
  palette = c("#E7B800", "#2E9FDF"),  # 双色区分突变/野生
  legend.title = "TP53 Status",
  legend.labs = c("Mutated", "Wildtype"),
  break.time.by = 5      # 规范时间轴刻度
)